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INTRODUCTION 

1 

UNC-SAM-2,  the  second  United  Nuclear  Corporation  Stochastic  Approximation 
Method,  is  a  code  system  in  the  FORTRAN  language  for  evaluating  the  time  de¬ 
pendence  of  the  transport  of  neutrons  or  photons  through  matter  using  Monte  Car¬ 
lo  techniques.'  '  N  V  „  '  ^  ‘  ' 

This  program  represents  the  latest  result  of  a  continuing  effort  to  develop  Monte 
Carlo  programs  which  will  handle  complex  geometrical  configurations  with  im¬ 
proved  versatility  and  speed,  -tt^is  based  upon  its  predecessor  codes  UNC  -SAM 1 
and  ADONIS.2  We  are  indebted  toMTfi.  Kalos  for  the  development  of  many  of  the 
numerical  techniques  employed  in  these  predecessor  programs  and  for  contribu¬ 
tions  to  some  of  the  improvements  which  have  been  incorporated  in  UNC-SAM-2, 
and  to  W.  Guber  for  major  contributions  to  the  formulation  of  UNC  SAM-2  as  well 
as  for  most  of  the  actual  coding  involved. 

Chapter  1  of  this  report  describes  some  of  the  new  importance  sampling  and  vari¬ 
ance  reduction  techniques  which  have  been  incorporated  into  the  code. 

Chapter  2  provides  a  short  description  of  the  UNC-SAM-2  program  and  each  of 
its  major  components. 

The  remaining  chapters  provide  more  detailed  descriptions  of  the  components 
of  the  code  with  representative  input  descriptions. 
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1.  METHODS  OF  ANALYSIS 


1.1  THE  BOLTZMANN  EQUATION 

1.1.1  General 

Let  us  define  the  density  of  particles  coming  out  of  collision  ip  (P,Q),  where  P 
represents  the  spatial  coordinates  X,  and  Q  the  direction  ft  and  the  energy  E. 

ip( P,Q)  satisfies  the  time -independent  Boltzmann  equation 

W?,0)  =  S(P,Q)  +  j>(P',Q')  T(P'-P  !Q')  C(Q'-QIP)  drp,  drQ,  (1) 

where  S(P,Q)  is  the  density  of  source  particles,  T(P'-PIQ')  is  the  transport 
kernel  which,  per  unit  particle  coming  out  of  collision  or  born  at  (P',Q'),  gives 
the  density  of  particles  coming  into  collision  at  (P,Q'). 

T  has  the  functional  form3 


T[X— X+Sfty(E,ft)]  =  ji(X  +  Sft',E)  e 


S 


M(X-+S'ft',E)dS' 


•  6{ft'*ftj)  6(ft'»ft2)  77(ft'-ft) 
where  ftt  and  ft2  form  an  orlhonormal  set  with  and 


(2) 


r?(X)  =  1  for  X  >  0 
i?(X)  =  0  for  X  <  0 

M(X, E)  =  the  macroscopic  total  cross  section  at  energy  E  of  the  material  in 
the  neighborhood  of  the  point  X. 

C(Q/‘^IP)  is  the  collision  kernel  which,  per  unit  particle  coming  into  collision 
at  (P  Q')>  gives  the  density  of  particles  coming  out  of  collision  at  (P,Q).  Details 
ol  , .  ble  forms  of  the  collision  kernel  handled  by  the  code  are  given  elsewhere. 
Here  we  only  assume  that  it  is  positive,  definite,  and  integrable  with  respect 
to  Q. 

The  flux  cp  (P,Q)  is  related  to  the  solution  ip(P,Q)  of  Eq.  1  through 

<p(P,Q)  =  fdrp,  ^(P/,Q)T(P,^PIQ)/m(P>Q)  (3) 

1.1.2  The  Estimates 

The  quantities  to  be  computed  are  of  the  form: 

H  =  Jv  <p(PfQ)  drp  drQ  (4) 

where  V  is  a  given  region  of  phase  space  or,  more  generally 

H  =  fv  <p(P,Q)  R(P,Q)  dtp  drQ  (5) 

where  R(P,Q)  is  a  response  function. 

0 

We  assume  here  that  R(P,Q)  is  definite  and  nonsingular,  that  the  integral  (Eq.  3) 
exists,  and  that  the  region  V  is  nondegenerate.  A  degeneration  of  a  specific 
type  [V(P,Q)  degenerating  to  a  single  point  P]  is  discussed  elsewhere. 

Substituting  Eq.  3  into  Eq.  4  we  can  write: 


2 


H  =  /  drp  fv  drp  drQ  ^Q)  TfP'-PIQVp^Q) 


(6) 


f 


1.1.3  Time  Dependence 


We  are  concerned  with  the  solution  of  the  time  independent  Boltzmann  equation. 
Even  so,  one  is  sometimes  interested  in  the  time  dependence  of  the  solution. 


The  source  Ji(P,Q,t)  is  assumed  to  be  a  known  function  of  time.  The  velocity  of 
particle  is  a  known  function  of  energy: 

v(E)  =  1.3833  x  106  VE  cm/sec  (E  in.  eV)for  neutrons 

v(E)  =  2.9978  x  1010  cm/sec  for  y  -rays. 


Collisions  are  assumed  to  take  place  instantaneously.  Under  these  assumptions 
the  Boltzmann  equation  becomes: 


4>(  P,Q,t)  =  S(P,Q,t)  + 


IP- PM 
v(E) 


T(P'-PIQ)  C(Q'—QIP)» 


•  dtp/  dTQ/ 


The  expression  of  the  flux  is: 


<p(  p 


,Q,t  )=f 


d rD/  4> 


P‘,Q,t  --^y  T(P'-PIQ)/m(P,Q) 


and  the  estimator  can  be  generalized  to: 


H  =  /y  <p(P,Q,t)  drp  drQ  dt 

where  V  is  now  a  region  of  phase  space  including  time. 

H  =  f  dTp./v  dTp  drQ  dt  4  P',Q,t  -  T(P'-P  IQ)/m(P,Q) 


1.2  THE  MODIFIED  BOLTZMANN  EQUATION 
1.2.1  General 


Let  us  introduce  an  importance  function  I  (P,Q),  or  a  weight  function  W(P,Q). 


I(P,Q) 


and  let  us  multiply  Eq.  1  by  this  importance  function 


ip( P,Q)  _  S(P,Q)  f  H P',0')  W(P,,Q/) 
W&Q)  WlP jQ)  +J  W(P',Q')  W(PQ') 


T(P'-PIQ') 


W(P,QQ 

WlP^Q) 


C(Q'-QIP)* 


(6a) 

*  dTjy  dTQ/ 


Let  us  define  a  modified  density  of  particles  coming  out  of  collision: 


4>  (P,Q) 


j/>(P,Q) 

w(p,q)  ’ 


(7) 


a  modified  source  density  function: 


S(P,Q) 


S(P,Q) 

wTp^q)  ’ 


(8) 


a  modified  transport  kernel: 

TtP'-PIQO  =  T(ly— PIQ  ), 


(9) 


and  a  modified  collision  kernel: 


C(Q'-QIP) 


_W(P,Q/) 

"  wrp,Q) 


C(Q/— QIP). 


(10) 


Substituting  these  quantities  into  Eq.  6a,  we  obtain  the  modified  Boltzmann  equa¬ 
tion: 

4 


4>( P,Q)  =  S(P,Q)  +  /^(P'.Q7)  T  (P'-P  IQ')  6  (Q'-Q  IP)  drp/  drQ/ 


(ID 


Substituting  them  into  Eq.  6,  we  obtain  a  new  expression  for  H: 

H  =  /  drp,  /v  dTpdrQ  4> (P',Q)  T  (P'-P  IQ)  W (P,Q)/p(P,Q)  (12) 

1.2.2  The  Importance  Function 

To  evaluate  the  desired  quantity,  H,  one  can  either  solve  Eq.  1  and  perform  the 
integration  (Eq.  5),  or  solve  Eq.  11  and  perform  the  integration  (Eq.  12). 

As  the  method  of  solution  is  the  Monte  Carlo  method,  the  estimate  of  H  will  be 
associated  with  a  certain  variance.  The  introduction  of  an  importance  (or  weight) 
function  gives  the  freedom  to  reduce  that  variance.  In  fact,  it  is  known  that,  if  the 
importance  function  is  taken  to  be  the  solution  of  the  problem  adjoint  to  Eqs.  1  and 
5,  the  variance  associated  with  the  estimate  (Eq.  12)  vanishes.  An  approximation 
to  this  “best”  importance  function  will  lead  to  finite  (but  small)  variance,  and 
should,  if  the  approximation  is  good,  lead  to  a  variance  smaller  than  the  one  as¬ 
sociated  with  the  unbiased  problem. 

In  the  code,  we  provide  the  possibility  to  take  into  account,  in  a  somewhat  crude 
manner,  the  spatial,  energy,  and  angular  dependence  of  the  weight  function. 

W  e  assume  that  the  weight  is  a  separable  function  of  energy  and  angle: 

w(x,e,3)  =wx(x)wE(x,E)wn(x,S5) 

where  the  first  term  is  not  really  necessary  but  is  introduced  because  of  com  - 
puter  coding  reasons. 

The  computer  program  is  restricted  to  assemblies  consisting  of  spatial  regions, 
each  containing  a  uniform  distribution  of  materials.  The  restrictions  for  the 
weight  function  that  the  code  can  handle  are  as  follows:  W-^fX)  is  constant  for 
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all  %  in  the  same  region;  Wg(X,E)  has  the  same  energy  dependence  for  all  X  in 
the  same  region.  The  energy  dependence  is  specified  as  follows.  A  set  of  energy 
bins  is  specified  by  the  energy  boundaries.  This  set  is  the  same  for  all  regions, 
but  can  berdifferent  from  any  other  energy  bins  (like  flux  bins).  For  each  spatial 
region,  a  set  of  numbers  is  specified,  corresponding  to  constant  energy  weight 
factors  in  each  energy  bin.  Similarly,  the  angular  weight  factor  W  ^(X,?2)  has 
the  same  angular  dependence  for  all  X  in  the  same  spatial  region.  The  angular 
weight  factor  is  assumed  to  depend  only  onto  =  where  %  is  specified  in 
each  region: 

Wfl(X,0)  =Wa,(X,S>fli)  . 


A  set  of  angular  bins  is  specified,  and  assumed  to  be  the  same  for  all  regions. 
For  each  spatial  region,  a  set  of  numbers  is  specified,  corresponding  to  constant 
angular  weight  factors  in  each  angular  bin.  There  are  three  different  ways  to 
specify  in  each  region.  The  simplest  is  to  specify  itself  as  a  constant 
vector  (aiming  vector)  or  to  specify  a  point  %  (either  target  or  source)  with 
ft*  specified,  at  each  point  X  in  the  given  region,  by 


X  -Rj 

llfc-Ril 


or  finally  to  specify  a  constant  angle  9 $  with  the  z-axis,  with  fij  specified,  at 
each  point  X  in  the  given  region,  by 


fiix  “ 


Vx2  +y2 


sin 


Oi..  =  — - -  sin  9i 

iy 


=  cos  0\. 
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At  present,  only  the  first  option  is  operational  in  the  code. 


1.3  SOLUTION  OF  THE  PROBLEM  BY  THE  MONTE  CARLO  METHOD 

1.3.1  General 

The  numerical  estimate  of  the  desired  quantity  H  is  obtained  as  follows:  Eq.  11 
is  solved  by  Monte  Carlo,  by  producing  a  population  of  points  (P,Q,t)  having  a 

a* 

density  proportional  to  ip  (P,Q,t).  The  estimation  of  H  is  then  carried  out  by 
performing  the  integration  (Eq.  12). 

The  solution  by  Monte  Carlo  goes  as  follows:  a  point  (P,Q,t)  is  generated  with 
a  probability  distribution  function  proportional  to  S(P,Q,t).  Points  where  par¬ 
ticles  enter  into  collision  are  generated  with  a  density  T  (Pold-^new  IQ)*  T^e 
time  of  flight  is  calculated.  Particles  getting  out  of  collision  are  generated  with 
a  density  C  (Qold*~Qnewlp)*  The  Process,  picking  from  T  and  then  from  C,  is 
repeated  until  no  further  possibility  exists  as  explained  below.  The  process  is 
then  repeated  for  another  source  point,  until  the  desired  number  of  source  par¬ 
ticles  has  been  exhausted. 

1.3.2  Picking  from  the  Modified  Source  Distribution  Function 

The  code  can  handle  either  a  source  tape  generated  by  an  external  source  genera¬ 
tor,  or  accept  as  input  a  description  of  the  source  under  some  restrictions,  and 
generate,  internally,  particles  from  that  source  modified  by  the  importance 
weights. 

The  source  tape  consists  of  a  collection  of  points  (P,Q,t)  with  a  specification  of  a 
region  number  in  which  they  occur,  and  a  set  of  two  weights:  Wp  and  Wc.  Wp  is 
the  statistical  weight  of  the  source  point.  Wc  is  a  normalization  factor  carried  by 
the  particle  and  by  the  entire  history,  with  all  scores  multiplied  by  that  factor. 

It  is  normally  set  as  one  (1).  Upon  picking  a  particle,  the  actual  weight,  W,  for 
the  particle  (P,Q)  is  computed  as  explained  in  Section  1.2.2,  and  the  factor 


F  =  Wp/W  is  calculated  and  interpreted  as  being  the  number  (which  may  be  frac¬ 
tional)  of  independent  source  particles  which  have  coordinates  (P,Q,t). 


The  built-in  source  generator  accepts  a  function  S(P,Q,t)  with  the  following  re¬ 
strictions: 

S(X,E,S2,t)  =  SX(X)  Se(E)  Sn(fi)  St(t). 

Sx(X)  is  assumed  to  be  constant  for  all  X  in  the  same  region;  power  densities 
must  be  specified  for  all  source  regions. 

Sg(E)  is  specified  by  a  table  of 


Sg(E)  dE  vs  E. 


Linear  interpolation  is  assumed  on 


S(E)  dE  vs  E. 


A  table  representing  the  Cranberg  fission  spectrum4  is  built  into  the  code,  but 
different  tables  can  be  read  in  as  an  option.  E high  ®low>  high  and  low  ener~ 

gy  cutoffs  of  the  problem  must  also  be  provided.  At  present,  Sfi(^)  can  be  either 
isotropic  or  monodirectional  only  (in  the  latter  case,  the  direction  of  the  source 
must,  be  specified).  St(t)  is  specified  by  a  table  of 

t 

St(t)  dt  vs  t. 


Linear  interpolation  is  assumed  between  entries  in  the  table. 
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Once  the  function  S(P,Q,t)  nas  been  specified,  the  code  will  pick  source  particles 
directly  from  a  probability  distribution  function  proportional  to  the  modified 
source  density  function  S. 


S(P,Q,t) 


S(P,Q,t) 
W(P,Q)  * 


The  mechanics  of  the  picking  routine  are  described  elsewhere.  Here  we  want 
only  to  mention  that  points  P,Q,t  are  generated  when  needed.  P,ft,  and  t  are 
picked  at  random  from  S,  whereas  the  energy  E  is  picked  in  a  stratified  way, 
in  descending  order  for  each  statistical  aggregate. 

As  we  pick  from  a  renormalized  modified  source  distribution,  answers  (fluxes) 
have  to  be  renormalized.  They  can  be  normalized  either  to  a  unit  source  par¬ 
ticle,  or  to  a  total  power.  The  normalization  is  achieved  by  multiplying  all  the 
region  weights  WX(X)  by  a  constant  normalization  factor,  once  and  for  all,  at 
the  beginning  of  the  problem. 

The  internal  picking  routine  delivers  a  point  (P,Q,t)  with  a  specification  of  the 
region  number  in  which  it  occurs,  and  a  set  of  two  weights:  Wp,  equal  to  the  ac¬ 
tual  weight  W  for  the  particle  (P,Q);  andWc  =  1.  The  factor  F  =Wp/W  is  then 
equal  to  unity,  and  the  particle  is  a  bona  fide  single  particle  with  coordinates 
(P,Q,t). 

1.3.3  Picking  from  the  Modified  Transport  Kernel  and  Scoring  of  Fluxes 

Once  a  particle  or,  more  generally,  F  particles  (P,Q,t)  are  picked  from  the  modi¬ 
fied  source  distribution  function  (or  from  a  modified  collision  kernel),  points 
where  particles  (P,Q)  enter  into  collision  must  be  generated  from  a  distribution 
which  is  F  times  the  modified  transport  kernel. 

Looking  at  Eq.  2  of  Section  1.1,  and  at  Eq.  9  of  Section  1.2,  the  points  P  are 
defined  by  P  =  P'  +  Sfi,  where  S  has  the  distribution  function  F(s): 
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F(s)  ds  = 


FW(g%Eg) 

W(X'+sS?,E?5) 


p$'+s^,E) 


expf 


p(X'+s'^,E)  ds' 


ds 


The  distribution  F(s)ds  is  not,  in  general,  a  probability  distribution  function,  as 
it  is  not  necessarily  normalized  to  unity.  The  original  transport  kernel,  as  well 
as 


are  normalized  to  unity.  But  the  biasing  function 


W(X',Efl) 

W(X'+sn,E^) 


in  general  destroys  that  normalization.  In  addition  the  factor  F  is  sometimes 
different  from  unity. 

Fig.  1(a)  shows  a  schematic  plot  of  the  function  F(s),  drawn  under  the  assumption 

that  points  into  a  direction  of  increasing  importance  (i.e.,  of  decreasing  weights). 

-  as 

As  can  be  seen,  the  s -dependence  is  e  in  each  region,  with  a  discontinuity  at 
each  boundary  where  a  change  in  weight  and  (or)  of  p  occurs. 

To  pick  s  from  such  a  distribution,  we  construct  the  function: 


G(s)  =  f  F(s)  ds 

•'O 


G(°°)  gives  the  normalization  of  F(s). 
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Fig.  1(a) 


If  N  is  an  integer  such  that  N  «  G(°°)  <  N  +  1,  we  will  arrange  to  pick  either  N 
or  N  +  1  distances  s,  with  probability  [N  +  l  -G(°o)]  or  [G^)  -N],  respectively. 
Furthermore,  the  collision  points  will  be  picked  in  a  stratified  way. 

To  pick  the  first  collision  point,  we  pick  a  random  number  ^  and  solve  [see 
Fig.  1(b)] 


G(sj)  = 

To  pick  the  second  collision  point,  we  pick  another  random  number  £2  anc*  solve 
G(s2)  =  1  +  £2>  etc* 

Before  solving  for  each  sn,  n  =  l,  ...,  n  +  l,  G(sn)  =  n-l  +  £n,  we  test  whether  such 
a  solution  exists,  i.e.,  whether  G(°°)  >  n-l  +  £n.  If  the  test  is  not  passed,  no  nth 
and  no  further  collisions  should  be  generated  from  that  flight. 

Actually,  we  do  not  pick  n  different  random  numbers,  but  just  a  single  one  at  the 
beginning  of  the  flight,  and  set  £2  -  £3  =  =  £n  = 

In  the  actual  coding,  the  tracking,  integration,  and  selection  of  collision  points  is 
done  simultaneously. 

In  describing  the  process,  it  is  implied  that  we  always  have  to  track  to  00  (or  to 
“escape”)  in  order  to  distribute  collision  points.  To  avoid  unnecessary  tracking 
into  regions  of  decreasing  importance,  we  test,  at  each  boundary,  the  integrand 
F(s)  vs  a  fixed  cutoff  value  F0  (a  number  less  than  1,  which  we  like  to  set  as  0.05). 
If,  at  some  s0,  F(s0)  becomes  smaller  than  F0,  we  Russian-roulette  the  remaining 
flight;  with  probability  [1-F(s0)]  we  cut  off  the  tracking,  whereas  with  probability 
F(s0)  we  multiply  the  remaining  density  function  F(s)  by  1/F(s0)  and  continue 
tracking.  The  case  is  illustrated  in  Figs.  1(c)  and  1(d). 
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The  time,  t,  at  which  the  interaction  is  entered  at  distance,  s,  is  given  by 
t  =  t'  +  s/v(E). 


Fluxes  are  being  scored  simultaneously  with  tracking.  As  seen  from  Eq.  12  of 
Section  1.2,  we  have  to  perform  the  integral 

H  =  /  drp,  /v  drp  drQ  i  (P',Q)  T  (P'-P  IQ)  W(P,Qy.u(P,Q) 

or,  if  we  are  interested  in  time  dependence,  the  integral 


/  dt'  fy  dt  drpdTQ  #(P',Q,t')  cJm') 


!P-P'I 

‘‘vW 


T(P'-PIQ) 


W  (P,Q) 

liTp^y 


Let  us  assume  that 


V  =  VpVEVt 


where  Vp  =  a  certain  region  of  space 
VE  =  a  certain  energy  interval 
Vt  =  a  certain  time  interval, 

the  expression  of  H  becomes 

H  =  f  dTpJdt'/dTjj/^dE  i  iJ(P',E,s5,t')/Vtdt/Vpdrp6[(t-t')  • 


FT(P'-P  IE,n)  • 


The  first  four  integrations  are  performed  by  Monte  Carlo  in  the  process  of  pick¬ 
ing  from  the  source  distribution  and  collision  kernels.  The  last  two  integrations 
are  performed  analytically. 
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dTp  dfft-O-i^gpJ  FT(P'-PIE,S)W(P,E,Syn(P,E)  = 


-fv t  dt/vs  ds  fi|(t-t')  -  |  F(s)  W(3^+s£,E,$)/p$'+s$,E) 


as  they  involve  the  same  integrations  needed  to  obtain  the  cumulative  distribu¬ 
tion  G(s). 

1.3.4  Picking  from  the  Modified  Collision  Kernel 

Given  a  particle  (P,Q',t)  going  into  collision,  we  are  supposed  to  pick  particles 
(P,Q,t)  from  the  modified  collision  kernel: 


C(<5'-QIP)  -  C(Q'— Q IP) 

Picking  is  not  done,  as  it  should  be,  directly  from  C,  but  rather  from  the  original 
collision  kernel.  However,  the  weight  W(  P,Q')  is  attached  to  the  particle.  Upon 
the  particle  being  picked  up,  the  fraction 

F  _W(P,Q') 
r  "  W(P,Q) 


is  calculated,  and  interpreted  as  being  the  number  (which  may  be  fractional)  of 
independent  particles  coming  out  of  collision  with  coordinates  (P,Q,t).  These 
will  be  split  or  Russian -rouletted  during  transport  as  we  have  seen  in  the  pre¬ 
ceding  section. 


The  possible  forms  of  the  collision  kernel,  and  methods  to  pick  from  it,  are  de¬ 
scribed  in  UNC-5093.1  The  present  code  uses  the  same  methods. 
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2.  SHORT  DESCRIPTION  OF  THE  CODE 


The  code  system  is  written  in  FORTRAN  for  the  Control  Data  Corporation’s  (CDC) 
1604-A  computer  and  for  the  Ballistic  Research  Laboratory’s  BRLESC  computer. 
A  32-K  core  and  from  2  to  8  magnetic  tapes  are  required.  UNC -SAM-2  is,  in 
reality,  a  chained  series  of  independent  programs  which  process  cross-section 
and  geometry  data,  do  the  transport  problem,  and  edit  the  results.  A  brief  sum¬ 
mary  of  the  independent  programs  contained  in  UNC -SAM-2  follows.  Fig.  2  indi¬ 
cates  the  information  flow  in  the  system. 

GENDA 

The  Generation  of  Data  routine  (GENDA)  starts  with  a  fundamental  set  of  neutron 
or  photon  cross-section  data  tabulated  as  a  function  of  energy.  It  allows  for  the 
basic  tabulation  of  the  total  cross  section,  the  scattering  cross  section,  the 
Legendre  expansion  coefficients  of  the  differential  scattering  cross  section,  the 
inelastic  scattering  cross  section  for  both  level  and  continuum  scattering,  etc. 
These  data  are  tabulated  generally  on  differing  energy  meshes.  GENDA  proc¬ 
esses  the  data  to  generate  a  set  of  cross-section  data  on  a  desired  cross-sec¬ 
tion  set  by  either  linear  or  logarithmic  interpolation.  The  specti  am  of  neutrons 
emitted  by  continuum  inelastic  scattering,  according  to  statistical  theory,  is  used 
for  heavier  materials.  For  lighter  materials,  such  as  Li  and  Be,  special  treat¬ 
ments  are  available.  The  output  of  this  routine  is  a  tape  which  may  be  used  as 
input  to  GENPRO. 
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Fig.  2  —  SAM -2  Flow  Diagram 
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GENPRO 


The  GENDA  Processor  (GENPRO)  routine  uses  '  he  GENDA  output  to  determine 
the  tables  of  probabilities  peculiar  to  a  nuclide,  required  for  a  Monte  Carlo  code, 
at  a  prescribed  energy  mesh.  For  instance,  the  probability  that  a  neutron  suf¬ 
fers  an  elastic  scattering  is  given  by  the  ratio  of  elastic-to-total-scattering  at 
that  energy.  The  probability  of  inelastic  scattering  from  one  energy  to  another 
is  computed  in  this  routine.  At  its  conclusion,  one  has  a  library  tape,  called  the 
Element  Data  Tape  (EDT),  on  a  given  energy  mesh  to  be  used  in  subsequent  proc¬ 
essing. 

Up  to  this  point,  no  problem -oriented  input  is  needed.  The  functioning  of  later 
routines  (TUNC,  GASP)  will  be  dependent  on  such  problem -oriented  data. 

Both  GENDA  and  GENPRO  codes  are  described  in  detail  in  UNC-5093.1 

TUNC 

TUNC  is  the  main  program  of  an  overlay  tape.  Fig.  3  indicates  the  information 
flow  in  the  system.  It  essentially  calls  in  sequence  the  following  subprograms. 

BAND 

Given  the  concentrations  of  the  materials  contained  in  the  configuration, 
an  Organized  Data  Tape  (ODT),  bearing  the  total  macroscopic  cross  sec¬ 
tions  as  well  as  other  probabilities  for  mixtures  of  nuclides,  is  created  by 
the  BAND  routine.  These  are  the  problem -dependent  cross-section  data 
tabulated  against  energy  at  the  prescribed  energy  mesh  previously  used  in 
GENDA  and  GENPRO. 

BEDIT 

Edits  an  existing  Organized  Data  Tape. 
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Fig.  3  —  TUNC  Flow  Diagram 
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GEOM 


This  is  a  routine  which  takes  a  simplified  geometrical  description  of  the 
physical  system,  as  provided  by  the  problem  originator,  and  produces  the 
rather  complex  set  of  data  required  by  the  transport  program. 

Three  types  of  geometries  are  available:  Type  1  is  for  axially  symmetrical 
problems;  Type  2  is  for  spherically  symmetrical  problems;  and  Type  3  is 
for  more  complicated  geometries,  in  which  case  the  problem  originator 
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must  decompose  space  into  a  set  of  boxes  known  as  “ordinary  regions.” 
Inside  these  boxes  may  be  placed  “nonordinary  regions”  consisting  of 
spheres,  cylinders,  wedges,  or  other  boxes.  Nonordinary  regions  may  en¬ 
close  other  nonordinary  regions. 

INPUTD 

Reads  in  the  remaining  problem  specifications,  namely,  the  assignment  of 
compositions  to  geometrical  regions,  the  assignment  of  regions,  energy, 
and  angular  dependent  statistical  weights,  and  the  assignment  of  scoring 
regions  (i.e.,  regions  where  average  fluxes  are  to  be  calculated),  as  well 
as  the  positions  of  point  detectors. 

The  source  specification  is  also  read  in.  The  source  can  be  defined  by 
either  an  external  source  tape,  or  by  a  description  of  the  energy  spectrum 
and  by  a  power  density  in  different  regions. 

MONTE 

Given  the  output  of  the  previous  subprograms,  one  is  now  ready  to  proceed 
with  the  transport  program  MONTE.  The  program  picks  source  particles, 
tracks  them  to  escape  (scoring  flux  contributions  where  needed),  distributing 
(and  storirg)  collision  points.  Particles  coming  out  of  collision  (if  any)  are 
subject  to  the  same  treatment  as  source  particles.  The  calculation  is  done 


in  statistical  groups  so  that  variances  can  be  calculated.  The  energy  range 
is  broken  up  into  supergroups  to  decrease  memory  limitations.  Infoima- 
tion  (cross  sections  and  scores)  relevant  to  only  one  supergroup  is  in  the 
computer  memory  at  any  time.  Answers  for  each  statistical  group  and 
each  supergroup  are  written  on  a  scratch  tape. 

PEDIT 

Reads  and  edits  the  output  tape  of  MONTE,  giving  fluxes  and  standard  devia¬ 
tions. 

TRANSMIT 

TRANSMIT  accepts  as  input  an  interaction  and  transmission  tape  generated  by 
MONTE.  It  selects  only  transmitted  particles,  i.e.,  particles  which  enter  a  given 
region  during  the  MONTE  tracking,  and  writes  out  a  source  tape  consisting  of 
these  particles.  It  also  edits  that  tape. 

GASP 

The  routine  Gamma  Source  Particle  generator  (GASP)  computes  gamma-ray  pro¬ 
duction  from  a  primary  neutron  problem  to  serve  as  a  source  for  a  secondary 
gamma-ray  problem.  GASP  requires: 

1.  An  interaction  tape  (generated  by  a  previous  TUNC  problem)  containing 
a  record  of  all  interactions  able  to  cause  gamma -ray  emission 

2.  Probability  numbers  for  each  nuclide  of  the  problem  giving  the  probable 
number  of  gamma  rays  emitted  as  a  function  of  energy. 

The  routine  can  cope  with  the  generation  of  other  secondaries  upon  modification 
of  the  probability  number  input.  The  output  of  either  GASP  or  TRANSMIT  is  a 
source  particle  tape  which  becomes  input  to  TUNC. 
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RESPONSE 


This  program  accepts  as  input  the  output  tape  of  MONTE.  It  calculates  integrals 
/  R(E)  (p  (E)dE  and  associated  statistical  errors  for  specified  functions  R(E)  in 
specified  regions. 


21 


3.  DATA  GENERATOR  AND  PROCESSOR  ROUTINES  (GENDA  AND 
GENPRO)  AND  THE  MAIN  CODE,  TUNC 

The  GENDA  and  GENPRO  routines  are  described  in  detail  in  UNC-5093.1  These 
codes  have  not  been  modified.  Nuclides  are  identified  on  input  to  GENDA  by  a 
10 -digit  integer  and  by  an  atomic  mass.  In  our  old  Monte  Carlo  codes,  the  only 
relevant  identification  of  the  element  was  the  integral  part  of  the  atomic  mass. 

In  the  present  code,  however,  the  relevant  identification  is  the  integer  consisting 
of  the  five  right-most  digits  of  the  10-digit  number  read  in  by  GENDA.  A  con¬ 
venient  way  to  identify  an  element  is  ZZAAA  where  ZZ  is  the  atomic  number  and 
AAA  the  mass  of  the  isotope  (set  to  zero  for  the  natural  element).  It  is  mentioned 
elsewhere  in  the  UNC-5093  report1  that  (n,2n)  reactions  are  ignored  in  the  inter¬ 
action  routines,  except  for  beryllium  for  which  two  neutrons  are  emitted  follow¬ 
ing  a  reaction  identified  as  inelastic.  This  is  still  true  in  the  present  code,  be¬ 
ryllium  being  recognized  as  having  an  identifier  04AAA  (any  AAA). 

TUNC  is  the  main  program  of  an  overlay  tape.  The  input  consists  of  two  cards, 
a  title  card  and  a  data  card  containing  IBAND  and  NBAND. 

If  IBAND  =  0  TUNC  will  call  BAND  which  creates  an  Organized  Data 

Tape  (ODT),  and  then  will  call  BEDIT  which  edits  the  ODT 

If  IBAND  =  1  TUNC  will  call  BEDIT,  which  edits  an  existing  ODT 

If  IBAND  =  2  TUNC  bypasses  both  BAND  and  BEDIT,  r.id  programs 

GEOM,  INPUTD,  MONTE,  and  PEDIT  are  called  sequen¬ 
tially. 
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4.  BAND 


The  purpose  of  BAND  is  to  arrange  the  cross-section  data  in  an  organized  fash¬ 
ion  for  the  particular  problem  in  process. 

The  GENDA  and  GENPRO  programs  are  used  to  generate  an  Element  Data  Tape 
(EDT)  which  contains  microscopic  cross  sections  and  probability  tables  for  all 
nuclides  of  interest.  This  EDT  can  be  considered  as  a  library  tape  for  any  sub¬ 
sequent  Monte  Carlo  problems. 

For  a  specific  problem,  however,  the  system  is  divided  into  regions  which  con¬ 
tain  given  nuclide  compositions.  Each  composition  may  be  described  by  a  set 
of  nuclides  and  their  respective  concentrations.  When  a  particle  enters  colli¬ 
sion  at  some  point  within  a  region  of  given  composition,  the  Monte  Carlo  program 
has  to  select  the  nuclide  in  the  composition  with  which  the  particle  collides. 


The  probability  that  a  particle  with  energy  E  will  interact  with  nuclide  k  of  the 
composition  R  is  given  by 


P(k,R) 


Ck,RaT^k»E^ 

^T/R^ 


where  Mt,r(e)  is  the  tota*  macroscopic  cross  section  of  the  composition,  and 
C^r  is  the  concentration  of  nuclide  k  in  region  R. 
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A  BAND  input  form  is  shown  in  Chapter  13.  The  BAND  output  is  a  long  array 
in  variable  dimension,  where  floating  cross-section  numbers  are  packed  together 
with  fixed  number  indices  in  single  words.  Chapter  5  describes  the  edit  of  the 
ODT. 
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5.  BEDIT 


BEDIT  edits  the  ODT  generated  by  BAND.  The  edit  is  by  energy  band. 
First  Line 

BAND-FLAG,  sequential  number  of  the  band  being  edited. 


Second  Line 

NO.  of  BANDS,  total  number  of  bands  in  problem. 

Third  Line 

Is  a  heading  for  the  table  that  follows. 

JEN1  is  the  location  of  the  first  word  of  the  energy  mesh  table; 

JEN2  is  the  location  of  the  next  to  the  last  word  of  the  energy  mesh  table. 

Next  Heading 

REGION  Composition  number 

LOCMU  The  location  of  the  first  word  of  the  macroscopic  total 
cross-section  table  for  that  composition 

LOCCON-  The  location  of  the  first  word  of  the  table  giving  the  con¬ 
centrations  of  different  nuclides. 

This  is  followed  by  a  table  for  all  compositions  present  in  the  problem. 
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Next  Heading 

CONC  Concentration  of  nuclide 

LOC  Location  of  first  word  of  microscopic  cross-section  data 
fox  the  nuclide 

ATWT  Atomic  weight  of  the  nuclide. 

This  is  followed  by  a  table  for  all  nuclides  in  each  composition  present  in  the 
problem. 

This  is  followed  by  a  printout  for  each  nuclide  present  in  the  problem. 

Heading 

LOELEV  Location  of  first  word  of  table  giving  discrete  energy 
levels 

ATWT  Atomic  weight  of  nuclide 

LPLEV  Number  of  discrete  energy  levels 

LCHI  Length  of  CHI -table 

LENN  Length  of  ENN-table 

The  meaning  of  CHI  and  ENN  tables  is  described  in  the  GENPRO  section  of  UNC- 
5093.  We  will  mention  here  that  the  angular  distribution  (E,  cos  9)  is  divided 
into  N  equally  probable  intervals  of  cos  9.  The  N+l  boundaries  of  cos  9  define 
N+l  values  of  x  =  (1  +  cos  9)/ 2,  which  are  stored  in  a  CHI  table  for  the  energy  E. 
The  length  of  the  CHI-table  is  N+l. 

Similarly  the  spectrum  of  inelastic  neutrons  P(E,E')  is  divided  into  M  equally 
probable  intervals  of  fi',  and  the  M+l  boundaries  of  E'  are  stored  in  an  ENN- 
table  for  the  energy  E.  M+l  is  the  length  of  the  ENN  table. 

Next  Heading 

E  Total  cross  section,  in  barns 

SIGMA  If  CHI  =  i,  elastic  scattering  is  isotropic  in  the  CM  system. 

If  =  2,  elastic  scattering  is  with  hydrogen  (special  routine). 
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SIGMA  If  =  4,  elastic  scattering  isotropic  in  the  lab  system. 

(Continued)  If  =  5,  the  Klein-Nishina  formula  is  to  be  used  (y- 

rays  only).  If  >  5,  then  CHI  is  the  location  of  the 
first  word  of  the  appropriate  CHI-table. 

SCATTER  Elastic  scattering  cross  section,  in  barns 

PLEV  If  =  0,  inelastic  scattering  is  not  via  discrete  levels. 

If  >  0,  PLEV  is  the  location  of  the  first  word  of  a 
table  giving  probabilities  to  excite  the  various  levels. 

ABSORB  Absorption  cross  section,  in  barns 

ENN  If  =  0,  inelastic  scattering  does  not  give  rise  to  a  con¬ 

tinuous  spectrum  of  inelastic  neutrons.  If  >  0,  ENN 
is  the  location  of  the  first  word  of  the  appropriate 
ENN  table 

L  The  location  of  the  total  cross-section  entry. 

This  is  followed  by  a  listing  of  all  the  tables  referred  to  (CHI,  ENN,  PLEV,  ELEV), 
together  with  their  relative  locations. 

After  the  printout  for  all  the  elements  comes  the  printout  of  macroscopic  cross 
sections  for  all  the  compositions,  in  cm”1  vs  E,  in  eV. 

This  is  the  end  of  the  edit  for  neutron  problems.  For  y-rays,  one  has,  in  addition, 
a  printout  of  EjCiZj  for  each  composition,  where  Cj  is  a  concentration  and  Z[  is 
an  atomic  number. 
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6.  GEOM 

Monte  Carlo  programs  track  simulated  particles  through  a  specified  geometrical 
configuration,  undergoing  all  the  interactions  that  an  actual  particle  is  expected 
to  undergo.  The  particle  flight  may  pass  through  media  having  different  proper¬ 
ties.  Therefore,  it  is  of  primary  importance  to  know  at  all  times  (1)  in  which 
medium  the  particle  is,  (2)  when  the  particle  leaves  the  medium,  and  (3)  which 
medium  it  will  be  entering.  Different  media  may  be  described  as  containing  given 
material  compositions  in  three-dimensional  geometrical  figures. 

However,  it  is  necessary  to  decompose  these  figures  into  particular  geometrical 
forms  which  can  be  handled  by  the  MONTE  program. 

The  GEOM  program  processes  the  input  which  consists  of  the  description  of  the 
desired  configuration  by  means  of  these  elementary  geometrical  forms  and  stores 
the  information  into  tables  for  rapid  data  access  in  subsequent  Monte  Carlo  pro¬ 
grams. 

Three  distinct  types  of  geometries  can  be  handled: 

1.  Axially  symmetric  geometries  (MASQUE) 

2.  Spherically  symmetric  geometries 

3.  Other  geometries  (UNC-SAM  geometry) 

The  input  common  to  all  three  types  of  geometries  is  a  title  card  (80  characters) 
followed  by  a  data  card  giving: 
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ITYPE  Type  of  geometry 

1  -  MASQUE 

2  -  SPHERICAL 

3  -  UNC-SAM 

NRGN  Number  of  distinct  regions 
NSRF  Number  of  distinct  surfaces 

IREX  Oi  ainal  number  of  the  region  declared  to  be  the  escape 
region  (>NRGN). 

6.1  AXIALLY  SYMMETRIC  GEOMETRY  (MASQUE) 

The  axis  of  symmetry  is  defined  to  be  the  z-axis.  The  geometry  is  specified 
in  the  (z,  positive  r)  half-plane.  Regions  in  this  half-plane  are  restricted  to  con¬ 
vex  quadrangles.  The  geometry  is  specified  by  first  defining  all  the  NSRF  bound¬ 
ing  lines,  and  then  defining  each  of  the  NRGN  regions  by  specifying  the  four  lines 
which  bound  them,  and  by  giving  information  on  the  neighboring  regions,  as  ex¬ 
plained  below. 

An  input  form  for  MASQUE  is  shown  in  Chapter  13. 

The  description  of  the  bounding  lines  is  as  follows.  First  one  numbers  all  the 
NSRF  lines  (including  the  z-axis),  from  one  (1)  to  NSRF,  in  arbitrary  order.  One 
then  defines  each  of  the  lines  in  sequential  order,  by  any  two  points  (R^Zj)  (R2,Z2) 
on  the  line.  The  input  for  each  line  consists  of  a  single  card  with  R1,Z1,R2,Z2,I, 
where  I  is  the  ordinal  number  of  the  line. 

All  the  quadrangular  regions  must  also  be  numbered  from  one  (1)  to  NRGN  in 
arbitrary  order.  One  then  has  to  define  each  of  the  regions,  in  sequential  order. 
For  each  region  IR,  a  card  gives 


L^L^Lj,  L4,J1,J2,J3,Ji,IR 


where  1^,  i  =  1,  2,  3,  4  are  the  ordinal  numbers  of  the  four  bounding  lines  of  re¬ 
gion  IR,  given  in  either  clockwise  or  counterclockwise  order.  If  there  is  a  single 
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neighbor  of  region  IR  across  line  Li,  then  Ji  is  set  equal  to  the  number  IR'  of 
that  neighboring  region.  If  there  is  more  than  one  neighbor,  say  nj  neighbors 
across  boundary  Li,  the  boundary  Li  is  called  a  complex  line  and  Ji  is  set  equal 
to  (—Hi).  If  the  region  IR  has  no  complex  lines,  the  input  for  that  region  is  com¬ 
plete.  If  complex  lines  exist,  further  input  is  required  for  each  of  the  complex 
lines,  in  order  of  increasing  i. 

To  specify  the  ni  neighbors  of  region  IR  across  the  complex  line  Li,  one  has  to 
specify  the  n*  region  numbers  IRj  and  the  coordinates  of  the  points  along  line  Li 
where  neighbors  change.  Actually,  as  these  points  lie  on  a  line  which  is  already 
defined,  only  one  of  the  two  coordinates  is  sufficient.  The  criterion  for  the  choice 
of  the  coordinate  (r  or  z)  is  that,  if  the  absolute  value  of  the  slope  of  in  the 
z  vs  r  system  is  less  than  one  (1),  the  r  coordinates  are  to  be  used.  The  z  co¬ 
ordinates  are  to  be  used  in  the  other  case.  Let  us  call  Q  the  coordinate  used. 

The  input  consists  of  n^  pairs  of  numbers 

(IR^Qj)  (IR£,Q2)  (iRApV 

where  Qt  <  Q2  <  ....,  Qn,  thus  signifying  that  regions  IR  and  IRj  touch  each  other 
along  the  segment  of  line  I4  between  points  Qj  and  (Qe  is  understood  to  be 
the  low  Q  vertex  of  line  L^). 

If  a  region  IR  is  adjacent  to  the  escape  region,  the  corresponding  IR'  should  be 
set  equal  to  the  escape  region  number  IREX.  Regions  that  include  the  z-axis  in 
three  dimensions  are  bound  by  the  z-axis  (and  three  more  lines)  on  the  r-z  plane. 
The  region  adjacent  to  such  region  across  the  z-axis  is  itself,  the  corresponding 
IR'  is  therefore  equal  to  IR.  The  quadrilateral  bounding  a  region  is  allowed  to 
degenerate  to  a  triangle,  but  a  dummy  fourth  line,  passing  through  a  vertex  and 
not  crossing  the  triangle,  must  be  included  in  the  specification.  Any  number  can 
be  used  to  specify  the  region  neighboring  the  triangular  region  across  such  a 
dummy  line.  Note  also  that  the  whole  configuration  must  be  convex,  completely 
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surrounded  by  the  escape  region,  and  that  each  point  within  the  configuration  must 
belong  to  some  region  IR  <  NRGN,  IR  *  IREX. 

6.2  SPHERICALLY  SYMMETRIC  GEOMETRY 

In  this  case  regions  are  defined  as  spherical  shells.  The  input  consists  of  NRGN 
radii  in  increasing  order,  which  define  the  NRGN+1  regions;  region  1  extends 
from  the  center  to  the  first  radius.  (The  input  NSRF  and  IREX  to  GEOM  is  ir¬ 
relevant  if  ITYPE  =  2.)  An  input  form  in  shown  in  Chapter  13. 

6.3  UNC -SAM  GEOMETRY 

In  the  UNC -SAM  geometry,  all  space  is  divided  into  nonoverlapping  rectangular 
parallelepipeds  designated  ordinary  regions.  Spheres,  right-circular  cylinders, 
rectangular  parallelepiped  or  right  wedges  (nonordinary  region)  with  arbitrary 
orientation  can  be  wholly  contained  within  an  ordinary  region  or  another  nonor¬ 
dinary  region.  By  suitable  combinations  of  these  regions  it  is  possible  to  repre¬ 
sent  many  complicated  configurations. 

This  geometry  is  described  in  detail  in  UNC-5093.1  An  input  form  in  shown  in 
Chapter  13. 
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7.  INPUTD 

The  purpose  of  this  routine  is  to  read  and  store  in  organized  fashion  the  remain¬ 
ing  problem  specifications.  An  input  form  is  shown  in  Chapter  13. 

The  first  card  contains  the  limits  of  the  problem  as  shown  below: 

NSTART  The  ordinal  number  of  the  first  history  to  be  treated  (usu¬ 
ally  1) 

NSTOP  The  ordinal  number  of  the  last  history  to  be  treated 

NSTAT  Number  of  histories  per  aggregate  for  statistical  scoring 
of  flux 

NRMAX  Number  of  distinct  geometry  regions 

NG  For  a  neutron  problem  set  NG  =  0;  for  a  gamma  prob¬ 

lem,  NG  =  1 

NT  Number  of  output  time  bins  (if  NT  =  0,  time  dependence 

is  bypassed) 

NOUT  Number  of  output  energy  bins 

NUMSC  Number  of  flux  scoring  regions 

NRWL  Number  of  distinct  region  weights 

IREX  Ordinal  number  of  the  escape  region 

IRT  Ordinal  number  of  the  transmission  region  (coordinates 

of  particles  crossing  a  boundary  into  that  region  will 
be  written  on  magnetic  tape  as  ‘‘transmitted  particles”)* 

The  next  cards  deal  with  point  detectors: 
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NDET  Total  number  of  point -detectors 

NDFAP  Number  of  detectors  for  which  finite  variance  estimator 
is  to  be  used.  Scores  for  the  remaining  (NDET-NDFAP) 
detectors  will  be  by  statistical  estimation  (more  details 
will  be  given  in  Section  8.4). 


The  next  NDET  cards  give  the  coordinates  x,  y,  and  z  of  the  NDFAP,  and  (NDET- 
NDFAP)  detectors. 


The  next  card  gives  other  limits  of  the  problem: 

ECUT  Low-energy  limit.  Particles  degrading  below  E  =  ECUT 
will  be  scored  as  “degrades”  and  their  tracking  and 
scoring  is  terminated.  This  energy  cutoff  should  be  set 
higher  or  equal  to  the  low-energy  limit  of  cross-section 
input. 

ETH  Irrelevant  in  the  present  version  of  the  code. 

TCUT  Remote  time  cutoff:  tracking  and  scoring  of  particles 
aging  beyond  a  time  T  >  TCUT  is  terminated. 

FZ  Minimum  allowed  value  of  F  (see  Chapter  2).  (Preferred 

value  FZ  =  0.05.) 

EHIGH  High-energy  cutoff  of  the  problem.  It  should  be  set  small¬ 
er  or  equal  to  the  high-energy  limit  of  cross-section 
input. 

The  next  cards  give  the  (NOUT+2)  output  energy  bin  boundaries.  The  first  and 
the  last  energy  boundaries  must  be,  and  some  (or  ail)  of  the  intermediate  bound¬ 
aries  can  be,  flagged  by  minus  signs.  The  energy  ranges  defined  by  such  flagged 
boundaries  is  defined  as  output  supergroups,  the  significance  of  which  will  be 
explained  in  Section  8.1.  The  first  energy  boundary  must  be  >  EHIGH.  The  last 
one  must  be  <  ECUT. 


The  next  cards  give  the  (NT+1)  output  time  bin  boundaries.  The  first  time  bound¬ 
ary  must  be  >  TCUT.  The  last  one  must  be  equal  to  0. 
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The  next  cards  give  the  NRWL  region  weights.  These  NRWL  region  weights  are 
in  effect  numbered  from  one  (1)  to  NRWL,  and  later  on  will  be  referred  to  by 
their  ordinal  number. 


The  next  NRMAX  cards  define  the  parameters  of  each  of  the  NRMAX  geometrical 
regions,  the  first  card  applying  to  region  1,  the  second  to  region  2,  etc.  Each 
card  gives: 


ISC  If  ISC  >  0,  fluxes  in  this  geometrical  region  will  be  scored 
in  the  scoring  region  with  ordinal  number  ISC.  It  is  pos¬ 
sible  to  have  the  same  scoring  region  number  assigned  to 
several  different  geometrical  regions.  If  ISC  =  0,  fluxes 
will  not  be  scored. 

NREG  Ordinal  number  of  the  composition  of  materials  present  in 
the  geometrical  region  described  by  the  card 

IRW  Ordinal  number  of  the  region  weight  assigned  to  the  geomet¬ 
rical  region  described  by  this  card 

IEW  Ordinal  number  of  the  energy  weight  set  assigned  to  the  re¬ 
gion  (if  0  -  no  energy  importance  in  the  region) 

IAN  Ordinal  number  of  aiming  vector  assigned  to  the  region  (if 
0  -  no  angular  importance  in  the  region) 

IANG  Ordinal  number  of  angular  weight  set  assigned  to  that  re¬ 
gion. 

After  these  (NRMAX)  cards,  the  different  energy  weights  must  be  defined.  Sev¬ 
eral  sets  of  distinct  energy  weight  sets  can  be  given.  One  defines  a  single  set 
of  energy  bins,  and  each  energy  weight  set  is  defined  by  a  set  of  corresponding 
weights  assumed  to  be  constant  in  an  energy  bin. 


The  next  card  gives: 

NEWL  Number  of  energy  bins  for  energy -importance  weights 
NEW  Number  of  distinct  energy  weight  sets. 


If  these  two  are  nonzero,  the  next  cards  give  the  (NEWL+1)  energy  boundaries, 
and  the  (NEW)  weight  sets.  The  directional  importance  is  specified  next: 
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NAIML  Number  of  distinct  aiming  vectors 
NUMA.NL  Number  of  angular  bins  for  angular  importance 
NUMANG  Number  of  distinct  angular  bin  sets. 

If  those  three  numbers  are  nonzero,  the  next  cards  give  the  x,  y,  and  z  projec¬ 
tions  of  the  (NAIML)  aiming  vectors,  the  (NUMANL+1)  boundaries  of  cos  6  defining 
the  angular  bins  for  energy  weighting  (the  first  number  must  be  £1.  and  the  last 
must  be  and  finally  the  NUMANG  sets,  each  set  consisting  of  (NUMANL) 

weights,  with  each  weight  assumed  to  be  constant  in  the  corresponding  angular 
bin. 

The  use  of  region,  energy,  and  directional  weights  is  explained  both  in  Chapter  1 
and  Section  8.3.  But  to  clarify  the  input  form,  we  will  recall  that  the  weight  of 
a  particle  at  a  given  X  flying  with  energy  E  in  direction  ft  is  given  by 

W(X,  E,  S?)  =  W£  (X)  We(E)  Wj5  (ft). 

Knowing  X  we  know  the  geometrical  region  number,  and  therefore,  IRW,  IEW, 

IAM,  and  JANG. 

IRW  is  the  ordinal  number  of  the  region  weight,  therefore,  (X)  is  known.  If 
IEW  =  0,  Wg  (E)  =  1;  if  IEW  >  0,  one  has  to  find  the  energy  bin  where  E  belongs, 
and  look  up  the  corresponding  energy  weight  Wg  (E)  from  the  energy  weight  set 
number,  IEW.  If  IAM  =  0,  Wj$  (S?)  =  1;  if  IAM  >  0,  one  looks  up  the  IAM^  vector 
^IAM>  computes  cos  e  =  ft*fti,  locates  the  angular  bin  where  this  cos  e  belongs, 
and  looks  up  the  corresponding  angular  weight  W^  (ft)  from  the  angular  weight 
set  number,  IANG. 

The  program  INPUTD  reads  in,  prmts  back,  and  stores  (in  packed  form)  all  the 
information  mentioned  above,  and  lays  out  the  memory  to  provide  storage  area 
for  the  different  scores  expected  in  the  Monte  Carlo  part  of  the  calculation.  It 
then  calls  the  subroutine  SOUCAL  which  reads  in  and  processes  the  source  data. 
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Source  Data 


An  input  form  is  shown  in  Chapter  13  for  SOURCE. 


The  first  card  contains: 


NSR  Number  of  regions  where  the  source  extends.  If  this  num¬ 
ber  is  0,  an  external  source  tape  is  expected,  and  no  fur¬ 
ther  input  is  required. 

IFLAG  If  IFLAG  =  0,  a  built-in  Cranberg  fission  spectrum  will 

be  used  for  the  energy  spectrum.  If  IFLAG  >  0,  the  spec¬ 
trum  is  specified  later  in  input.  It  will  be  defined  by 
IFLAG  entries. 

ISW  Switch  determining  the  normalization  of  problem.  If  it  is 
0,  fluxes  will  be  normalized  to  a  single  (unbiased)  source 
particle.  If  it  is  1,  fluxes  will  be  normalized  to  the  “to¬ 
tal  power”  of  the  source  as  defined  below. 


This  is  followed  by  (NSR)  cards  specifying  the  (NSR)  source  regions.  Each  of 
these  cards  gives: 


ISR  Specification  of  a  geometrical  region  number 

P  Power  density  in  that  region 

ISO  Flag  indicating  angular  distribution: 

ISO  =  0,  isotropic 

ISO  =  1,  monodirectional.  If  ISO  =  1  in  any  of  the  (NSR) 
regioas,  a  monodirectional  source  will  be  assumed  in 
all  the  (NSR)  regions. 


After  this  follow  the  (IFLAG)  cards  specifying  the  energy  spectrum  assumed  to 
apply  to  all  the  source  regions.  This  is  a  table  of  E  vs  F(E)  where 


roo 

S(E)dE. 


Linear  interpolation  is  assumed  on  E  vs  log  [F(E)j.  The  first  entry  must  be  for 
E  >  Efoigh,  and  the  last  for  E  <  ECU£. 
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If  there  is  time  dependence  (NT  >  0),  the  time  dependence  of  the  source  must  be 
specified.  A  card  gives: 

NOT  Number  of  cards.  This  card  is  followed  by  NOT  cards  giving 
t  vs  f  St(t)  dt. 

Jo 

Linear  interpolation  is  assumed  between  entries  in  the  table. 

Finally,  if  the  source  is  monodirectional,  the  projections  £2X,  fty,  f2z  of  the  direc¬ 
tion  must  be  given  on  the  last  card. 


The  subroutine  SOUCAL  reads  in  all  this  input,  prints  it  back,  and  pre -computes 
tables  to  pick  directly  from  the  biased  source  distribution.  The  code  first  pre¬ 
computes  a  table  of 


SPEC  (I)  = 


S(E)dE, 


where  the  Ej’s  take  the  values  of  E^igh?  of  all  the  energy  boundaries  where  the 
energy  weight  changes,  and  Ecut.  The  cod<~  then  runs  through  all  the  source  re¬ 
gions,  and,  for  each  new  energy -importance  set  encountered,  pre-computes  a 
table  of 


SPEC  (I,  J) 


S(E)dE 

wJIeT 


where  J  runs  from  1  to  the  total  number  of  different  energy -importance  sets 
encountered  in  the  source  regions;  also  for  each  new  angular  importance  set  en¬ 
countered,  a  table 


P(I,K)  = 


dw 
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where  the  Wj’s  take  the  values  of  cos  6  at  which  the  angular  weight  changes,  and 
K  runs  from  1  to  the  total  number  of  different  angular  importance  sets  encoun¬ 
tered.  The  different  tables  are  renormalized,  and  both  the  modified  and  unmod¬ 
ified  integrated  source  are  computed  in  each  source  region.  The  former  quanti¬ 
ties  are  proportional  to  the  probability  with  which  particles  should  be  picked  in 
different  source  regions.  A  table  SOUR(L)  is  built  up,  which  gives  the  cumula¬ 
tive  probability  for  a  source  to  be  picked  in  the  region  for  M  >  L. 
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8.  MONTE 


8.1  THE  MAIN  PROGRAM 

The  program  MONTE  is  the  main  program  of  an  overlay  where  all  the  tracking 
and  collision  mechanics  are  performed.  The  program  arranges  the  c?dculations 
in  supergroups.  We  have  seen  that  the  program  BAND  arranges  cross  sections 
in  certain  energy  bands.  The  output  energy  bins  are  also  arranged  in  certain 
output  supergroups.  Cross-section  input  corresponding  to  a  single  band  can  be 
stored  in  computer  memory  at  any  given  time.  Scores  corresponding  to  a  single 
output  supergroup  can  be  stored  in  computer  memory  at  any  given  time.  The 
two  meshes  (bands  and  output  supergroups)  do  not  necessarily  coincide.  A  com¬ 
bined  mesh  defines  a  set  of  supergroups. 

The  program  Monte  starts  by  reading  in  the  highest  energy  cross-section  band, 
and  by  arranging  the  memory  for  the  highest  energy  output  supergroup.  The  high¬ 
est  of  the  low-energy  bounds  of  these  energy  ranges  defines  EBL,  the  low  energy 
of  the  supergroup  currently  treated. 

The  program  then  calls  the  source -picking  routine  SOUPIC,  and  examines  the 
energy  E.  If  E  <  EBL,  the  particle  is  stored  as  i  latent  by  calling  the  subroutine 
STORE,  and  SOUPIC  is  again  called.  When  E  >  EBL,  the  subroutine  CARLO  is 
called.  The  subroutine  CARLO  tracks  the  particle,  scores  contributions  to  fluxes 
when  needed,  and  distributes  collisions  along  their  track.  Particles  coming  out 
of  collR'on  are  also  tracked  if  their  energy  is  above  EBL;  particles  coming  out 
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of  collision  with  E  <  EBL  are  stored  as  latents  by  calling  STORE.  Control  is 
finally  returned  to  MONTE,  which  proceeds  to  the  next  source  particle,  until  a 
complete  statistical  aggregate  of  particles  has  been  treated  for  the  highest  super¬ 
group. 

At  this  point,  MONTE  switches  to  the  next  supergroup  by  either  reading  a  new 
band  of  cross-section  data,  or  by  writing  out  on  tape  the  set  of  scores  obtained 
and  preparing  the  memory  layout  for  the  next  supergroup,  or  both.  It  then  pro¬ 
ceeds  to  call  a  subroutine  PICK,  which  picks  latents  from  previous  supergroups. 

If  E  <  EBL  the  particle  is  stored  again  as  a  latent  by  calling  STORE.  If  E  >  EBL, 
CARLO  is  called.  The  procedure  continues  until  all  latents  have  been  examined, 
at  which  point  MONTE  switches  to  the  next  supergroup  until  the  low-energy  cut¬ 
off  is  encountered.  When  this  occurs,  MONTE  switches  to  the  highest  supergroup, 
and  proceeds  to  treat  the  next  statistical  aggregate  of  particles.  The  calculation 
terminates  when  a  history  number  exceeds  the  cutoff  value  NHIST  specified  on 
input.  A  “blank”  interaction  record,  with  NHIST  =  NHIST  +  1,  is  written  on  the 
interaction  tape  and  all  tapes  are  rewound.  Control  is  transferred  to  TUNC  which 
calls  the  editing  program  PEDIT. 

8.2  SOUPIC  -  THE  SOURCE -PICKING  ROUTINE 

This  is  a  subroutine  which  picks  particles  from  the  biased  source  distribution. 

If  an  external  source  tape  is  to  be  used,  groups  of  35  source  particles  are  read 
from  tape  into  a  buffer,  and  returned  one  by  one  to  the  main  code.  The  quantities 
describing  a  source  particle  are: 

XE  ~  coordinates  of  the  particle 
IR  -  region  number  where  particle  is  born 
WB  -  direction  of  the  particle 
T  -  time  at  which  particle  is  born 
E  -  energy  of  the  particle 

NHIST  -  history  number  attached  to  the  source  particle 
F  -  statistical  weight  of  the  particle  generated 
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Wc  -  carry -along  statistical  weight  of  the  particle  (usually  set  equal 
to  unity). 

A  first  random  number  |  is  compared  to  the  table  SOUR(L).  The  smallest  L  for 
which  SOUR(L)  >  £  determines  the  region  IR  =  ISR(L)  to  be  picked  from.  Stand¬ 
ard  techniques  are  used  to  pick  coordinates  of  points  uniformly  distributed  in  a 
region;  coding  is  included  for  any  region  in  either  spherical  or  MASQUE  geom¬ 
etry,  for  rectangular  parallelepipeds  (both  ordinary  and  nonordinary),  and  spheres, 
without  internal  bodies,  for  the  UNC-SAM  geometry. 

The  energy  is  picked  next.  A  stratified  random  number  (called  CE  in  the  code) 
is  obtained  (stratification  is  done  for  each  statistical  aggregate  of  source  parti- 

•r 

cles).  A  biased  random  number  £  is  then  obtained  by  interpolation  in  the  SPEC  vs 
SPEC(J)  tables  pre-computed  by  SOUCAL.  (The  J  is  determined  by  the  region 
number.)  Finally  the  energy  is  determined  by  solving  the  equation 


S(E)dE 


using  semi-log  interpolation  as  explained  in  Chapter  7. 
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The  direction  WB  of  the  source  particle  is  determined  as  follows.  If  the  source 
is  isotropic  and  there  is  angular  importance  sampling,,  the  cosine  of  the  angle 
between  the  particular  aiming  angle  and  the  direction  is  chosen  by  picking  a  ran¬ 
dom  number,  and  interpolating  between  the  angular  mesh  supplied  on  input  vs  the 
table  P(I,K)  pre-co  iputed  by  SOUCAL.  (The  K  is  determined  by  the  region  num¬ 
ber.)  A  random  azimuth  is  then  picked,  which  completes  the  specification  of  the 
direction.  In  the  absence  of  angular  importance  in  the  source  region,  standard 
techniques  are  used.  The  case  of  a  monodirectional  source  also  can  be  handled, 
provided  there  is  no  angular  importance  in  the  source  region.  Finally,  if  time 
dependence  is  to  be  determined,  a  time  T  is  determined  by  another  random  num¬ 
ber  £  and  the  solution  of  the  equation 

£  =  f  St(t)dt. 

•'0 

o  _  ^ 

The  quantities  communicated  to  the  main  code  are  XB,  IR,  WB,  T,  E,  NHIST,  F, 
and  Wc  (where  Wc  =  1),  and  F  is  the  weight  in  region  IR  at  energy  E  in  the  di¬ 
rection  WB. 

8.3  CARLO  -  THE  TRACKING  AND  SCORING  ROUTINE 

Source  particles,  either  read  from  source  tape  or  generated  internally  by  the 
subroutine  SOUPIC,  as  well  as  latents,  are  tested  on  energy.  If  the  energy  is  be¬ 
low  the  lower  bound  of  the  supergroup  currently  treated,  they  are  stored  as  la- 
tents.  For  energies  within  the  supergroup,  they  are  transmitted  to  the  subroutine 
CARLO  which  tracks,  distributes  further  collisions,  if  any,  and  scores  answers. 

Given  the  region  number  IR,  the  energy  E  and  the  direction  of  flight  WB,  a  total 
weight 


W  -  WX  x  WE  x  Wj5 
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is  calculated  and  a  new  F  given  by  the  ratio  F/W  is  treated  as  the  total  number 
of  biased  particles  born  simultaneously  at  XB,  T,  with  an  energy  E,  in  the  di¬ 
rection  WB.  The  carry -along  statistical  weight,  Wc,  can  be  considered  as  a  nor¬ 
malization  factor,  to  be  transmitted  to  all  future  collided  particles  due  to  the 
present  source  particle.  For  internally  generated  particles  the  quantity  F  will 
always  be  unity.  Externally  generated  particles  will  have  an  F  different  from 
unity  if  they  were  not  picked  from  the  properly  biased  source  distribution.  The 
quantity  F  is  tested  vs  an  input  cutoff  value  FZ(<1).  If  it  is  smaller,  a  game  of 
chance  is  played,  and,  with  probability  (1-F)  the  particle  is  killed;  whereas  with 
probability  F  the  particle  is  kept  with  F  set  as  unity.  These  two  events  are 
scored  as  either  (F  x  W)  “kills”  or  (l-F)xw  “births.” 

If  the  particle  survives  the  test,  a  subroutine  DR  is  called,  which  provides  the 
total  macroscopic  cross  section  p  in  the  region  IR  (this  routine  is  almost  identi¬ 
cal  to  the  DR  routine  described  in  UNC-5093).1  A  random  number  ^  is  picked 
for  picking  future  collision  points. 

A  geometry  routine,  Gl(Si, IR',X'),  is  call'-  to  provide  the  distance  St  to  the 
first  region  boundary  encountered  from  XB  in  the  direction  WB.  IR'  is  the  re¬ 
gion  number  on  the  other  side  of  the  boundary,  and  X'  is  the  point  of  intersection 
of  the  track  with  the  boundary.  If  the  region  IR  is  a  scoring  region,  the  contri¬ 
bution  to  the  flux  is  calculated  and  scored,  as  explained  in  Chapter  1.  Points 
where  particles  enter  into  collision  are  picked  by  computing 

{'  =  I,  -  F(l-e"MSl). 

If  <  0,  a  particle  comes  into  collision  in  region  IR  at  X  =  XB  +  S :ti  where  S'  = 
-1/p  log  (1  -  li/F);  the  interaction  routine  DR  is  called  to  pick  a  particle  coming 
out  of  collision.  The  interaction  is  written  on  tape,  and  the  particle  coming  out 
of  collision  (if  ever)  is  stored  as  a  latent.  is  then  increased  by  one,  and  con¬ 
trol  is  returned  to  the  part  of  the  code  which  computes  to  attempt  to  produce 
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more  collisions  in  region  IR.  When  >  0,  no  more  collisions  are  produced  in 
region  IR.  F  is  then  set  to  F  x  e  ^xw,  where  W  is  the  weight.  A  test  is 
made  whether  the  next  region  is  a  transmission  region  (if  it  is,  the  coordinates 
X',  IR,  energy,  time,  etc.  are  written  on  tape)  and  whether  it  is  the  escape  re¬ 
gion.  If  it  is  not  an  escape  region,  the  particle  tracked  is  moved  to  the  bound¬ 
ary  by  setting  XB  =  XP,  IR  =  IR',  etc.,  by  computing  the  new  weight  W,  and  by 
setting  a  new  F  equal  to  F/W.  If  F  <  FZ,  the  particle  is  either  killed  or  F  is 
set  equal  to  1,  both  events  being  properly  tallied.  For  F  <  FZ,  £  is  set  equal 
to  £',  the  DR  routine  is  called  to  obtain  the  new  total  cross  section  p,  and  con¬ 
trol  is  transferred  to  the  part  of  the  code  which  calls  the  geometry  routine 
G1  (S1?IR',X'),  thus  the  tracking  continues  until  either  a  kill  occurs  or  the  es¬ 
cape  region  is  reached.  In  the  latter  case  F  “escapes”  are  tallied. 

We  mentioned  above  that  particles  coming  out  of  collision  (if  ever)  were  stored 
as  latents.  An  absorbed  particle  does  not  come  out  of  collision,  but  is  scored 
as  W  “absorbs.”  Elastic  or  inelastic  particles  come  out  with  some  energy  E. 

If  E  is  less  than  the  low-energy  bound  of  the  supergroup  currently  treated,  the 
particle  is  stored  until  the  next  supergroups  are  treated  by  calling  the  subroutine 
STORE.  If  the  energy  is  high  enough,  the  particle  is  stored  in  a  local  buffer. 

Upon  termination  of  a  tracking,  the  subroutine  CARLO  examines  this  local  buffer, 
and,  if  any  latents  are  available,  picks  one  and  performs  the  tracking.  When  the 
buffer  is  empty,  control  is  returned  to  the  MONTE  program. 

8.4  FAP  -  THE  FLUX -AT-A -POINT  ROUTINE 

The  question  of  estimation  of  flux-at-a-point  by  Monte  Carlo  methods  was  inves¬ 
tigated  in  detail  by  M.  H.  Kalos.6 

The  methods  he  proposed  deal  with  obtaining  a  distribution  of  points  where  par¬ 
ticles  enter  into  collision  with  normal  Monte  Carlo  methods,  and  estimating  the 


once -more -collided  flux,  at  the  point  of  interest,  using  a  special  distribution  of 
the  last  collision  point. 


Assuming  a  collision  point  at  S,  a  detector  at  D  as  shown  in  the  sketch,  the  ex¬ 
pression  for  the  once-more-collided  flux  is 


where  gg(cu)  dw  is  the  probability  of  scattering  through  an  angle  whose  cosine  is 
in  the  range  (u),w+dco). 

For  isotropic  source  particles  one  has  to  estimate  in  addition  the  uncollided  flux. 

One  of  the  methods  Kalos  proposes  is  to  pick  the  intermediate  points  A  from  a 
distribution  given  by 
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(14) 


the  estimator  becoming 
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(15) 
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fl  =  16R  gS*w)  gA^cos  e 
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This  is  the  method  used  in  the  present  code.  It  turns  out  that  there  is  a  very 
simple  way  to  pick  intermediate  points  A  from  the  distribution  (Eq.  14),  if  a 
proper  choice  of  coordinate  system  is  made. 

Let  us  define  the  point  A  by  the  three  angles  f)u  02,  anc*  where  0t  and  02  are 
defined  on  the  sketch  and  cp  is  an  azimuth  around  SD. 

The  distribution  (Eq.  14)  becomes: 


R 

P 


9  x* 

r\  sin  6i  d0x  dcp  — *  d02 

- zrzz - - - 

r  1  r2 


Using  the  relationships: 


(16) 


ri  _  r2  .  R 
sin  (e2-  0i)  sin  0j  sin  02 

one  has 


r  =  sin  e>  n 
2  sin  e2 


and 


9r1  _  sin  R 
902  sin2  02 


Substituting  into  Eq.  16  one  obtains  the  distribution 
pr  d0j  d02  dcp 


(17) 


i.e.,  a  random  distribution  of  0\,  02,  cp  in  the  range 
0  <  0j  <  02  «  7T 


(18) 
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0  <  cp  <  2n 


(19) 


Once  6i  and  <p  are  picked,  the  cosine  of  the  angle  of  scattering  w  is  easily  ob¬ 
tainable,  and  the  score  (Eq.  15)  can  be  evaluated. 

In  the  actual  program,  coding  is  provided  for  either  the  procedure  outlined  above, 
which  gives  a  finite  variance  provided  the  detector  is  not  in  the  immediate  vicinity 
of  the  source  region,  or  for  a  simple  statistical  estimation  method,  which  gives 
finite  variance  if  the  detector  is  not  in  the  immediate  vicinity  of  either  the  source 
region  or  of  a  scattering  material.  In  that  case  the  quantity 

lids 

es  fcos  e)  - 

is  scored  for  each  collision  point  S,  where  6  is  the  angle  between  the  directions 
H  and  SD.  The  type  of  calculation  to  be  performed  is  specified  on  input,  by  spec¬ 
ifying  the  total  number  of  detectors  ND,  and  the  number  NDFAP  of  detectors  for 
which  the  first  procedure  is  to  be  applied.  The  remaining  (ND-NDFAP)  detectors 
will  be  treated  by  statistical  estimation. 

If  ND  *  0,  the  flux-at-a-point  routine  FAP  is  called  for  each  source  particle  by 
the  MONTE  program,  and  for  each  collision  point  by  the  CARLO  routine,  except 
if  inelastic  scattering  occurs  giving  an  outgoing  energy  below  the  supergroup 
being  currently  treated.  In  that  case  the  particle  is  stored  as  a  latent,  and  MONTE 
will  call  the  FAP  routine  with  that  latent  when  the  proper  supergroup  will  be 
treated. 

Once  called,  the  subroutine  FAP  examines  whether  the  particle  to  be  treated  is 
a  source  particle,  a  particle  coming  out  of  inelastic  scattering,  a  particle  coming 
out  of  elastic  scattering,  or  a  latent  from  a  previous  FAP  calculation.  In  the  first 
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three  cases,  the  trigonometric  functions  of  0t  and  e2  picked  according  to  the  dis¬ 
tribution  (Eq.  17)  are  obtained.  This  is  done  rather  efficiently  by  picking  directly 
the  sin  and  cos  functions  of  two  random  angles,  and  then  determining  which  cor¬ 
responds  to  the  smallest  angle,  and  therefore  which  should  be  assigned  to  0t;  the 
others  are  assigned  to  e2.  6\  and  02  are  kept  the  same  for  all  FAP -detectors  in 
the  problem. 

For  source  particles,  the  uncollided  flux  is  computed  (the  source  is  assumed  to 
be  isotropic). 

The  contribution  of  the  once-more-collided  flux  is  computed  as  follows.  The 
distances  rj  and  r2  are  given  by  r2  =  R  sin  flj/sin  02  and  rj  =  R  sin  -  r2  sin  e2. 

A  random  azimuth  around  SD  is  picked.  For  both  source  particles  and  inelastic 
particles,  gg  =  1,  and  the  outgoing  energy  is  as  set  by  the  main  code  (either  en¬ 
ergy  of  the  source  particle,  or  energy  of  the  inelastic  particle).  (No  energy - 
angular  correlation  is  assumed  in  the  code.  The  difference  between  lab  and  cen- 
ter-of-mass  system  is  also  ignored  for  inelastic  scattering.)  For  particles  com¬ 
ing  out  of  elastic  scattering,  the  angle  of  scattering  is  computed  in  the  lab  system. 
A  subroutine  GE  is  called.  For  neutrons,  it  performs  the  transformation  to  the 
CM  system,  looks  up  the  probability  gg,  and  computes  the  outgoing  energy  E.  For 
gamma  rays,  gg  and  E  are  given  by  the  Klein-Nishina  formula.  If  this  energy 
is  below  the  supergroup  currently  treated,  the  particle  is  stored  as  a  latent  of 
the  first  kind  jto  be  picked  up  and  returned  to  this  point  of  the  coding  at  a  later 
time  in  the  calculation. 

The  subroutine  DR  is  called,  and  the  total  cross  section  p  is  made  available.  A 
subroutine  TRALA  is  called,  which  tracks  the  particle  a  distance  rt  in  the  direc¬ 
tion  of  flight,  with  an  energy  E,  and  computes 

ffl 

Xt  =  I  pds; 

Jo 
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if  the  particle  escapes,  an  “infinite”  At  (At  =  27)  is  returned,  no  score  is  made, 
and  the  next  (if  any)  detector  is  treated.  If  the  particle  did  not  escape,  the  sub¬ 
routine  DR  is  called  to  determine  the  type  of  interaction  occurring  at  the  inter¬ 
mediate  point.  Absorption  corresponds  to  no  score.  For  inelastic  particles,  E  is 
picked  from  the  spectrum  by  the  DR  routine,  and  =  1.  For  elastic  scattering, 
the  subroutine  GE  determines  g^  and  E.  If  the  energy  E  is  below  the  supergroup 
currently  treated,  the  particle  is  stored  as  a  latent  of  the  second  kind,  to  be  picked 
up  and  returned  to  this  point  of  the  coding  at  a  later  time  in  the  calculation.  Fi¬ 
nally,  the  subroutine  TRALA  provides 

Jr*2 

I  pds, 

0 

and  the  total  contribution  is  calculated  and  scored,  and  the  next  detector  (if  any) 
is  considered,  unless  the  calculation  dealt  with  a  latent  of  either  first  or  second 
kind,  which  applies  only  to  a  particular  detector.  When  all  FAP  detectors  have 
been  considered,  the  statistical  estimation  detectors  are  treated  by  using  the 
coding  described  above  for  the  tracking  from  intermediate  point  to  the  detector, 
where  the  intermediate  point  is  replaced  by  the  original  collision  point. 

8.5  PICK  AND  STORE  -  THE  TREATMENT  OF  LA  TENTS 

We  have  seen  throughout  the  previous  sections  that  particles  degrading  below 
the  energy  EBL,  low-energy  limit  of  the  supergroup  currently  treated,  were 
stored  as  iatents  by  calling  the  subroutine  STORE.  They  were  later  picked  up 
by  MONTE  by  calling  the  subroutine  PICK.  The  two  subroutines  are  actually  a 
single  one  with  two  entry  points.  The  subroutine  INPUTD  allocates  the  memory 
to  data,  scores,  etc.  The  remaining  memory  is  assigned  to  the  subroutine  PICK- 
STORE,  to  be  used  as  a  buffer  for  Iatents.  One  “end”  of  the  buffer  is  assigned 
to  “degraded”  particles,  this  is  the  end  of  the  buffer  where  particles  are  being 
stored.  The  other  “end”  of  the  buffer  is  assigned  to  “unsorted”  particles,  i.e., 
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the  particles  to  be  picked.  Associated  to  each  end  of  the  buffer,  there  is  a  mag¬ 
netic  tape  to  be  used  when  the  buffer  overflows.  There  are  two  modes  of  opera¬ 
tion.  In  one  mode,  the  top  of  the  buffer  is  unsorted  and  the  bottom  is  sorted.  When 
the  switch  is  made  from  one  supergroup  to  the  next,  the  “unsorted”  part  is  empty, 
and  the  “degraded”  part  may  have  particles  which  become  “unsorted”  for  the 
supergroup  about  to  be  treated.  The  designation  of  the  buffers  (and  of  the  tapes) 
is  therefore  switched. 

There  is  no  set  boundary  between  the  two  “ends”  of  the  buffers.  The  number 
of  particles  in  the  “unsorted”  buffer  keeps  decreasing,  whereas  the  number  of 
particles  in  the  “degraded”  buffer  keeps  increasing,  and  can  increase  faster  than 
the  other  number  decreases.  Therefore,  the  two  parts  of  the  buffer  can  meet 
causing  an  overflow  of  the  buffer. 

It  is  then  determined  which  “end”  of  the  buffer  is  longest,  and  a  number  of  par¬ 
ticles  exactly  equal  to  one-half  the  total  length  of  the  buffer  are  written  from  the 
longest  “end”  onto  the  corresponding  magnetic  tape. 

When  the  “unsorted”  buffer  becomes  empty,  a  lest  made  whether  any  “un¬ 
sorted”  particles  are  available  on  the  corresponding  magnetic  tape.  If  none  are 
available,  the  calculation  has  been  completed  for  the  current  supergroup.  If  some 
are  available,  they  are  read  into  the  buffer  if  room  is  available.  If  room  is  not 
available,  it  is  made  available  by  writing  out  part  of  the  other  buffer  on  the  other 
tape;  the  length  of  the  record  written  out  from  one  “end”  is  equal  to  the  length 
of  the  record  to  be  read  into  the  other  “end.” 

The  subroutine  PICK-STORE  deals  with  different  kinds  of  latents.  The  quantities 
stored  are:  X,  E,  IR,  T,  I,  F,  NHIST,  WC,  and  J 12345,  where  X  is  the  position, 
ft  the  direction,  E  the  energy,  IR  the  region  number,  T  the  time,  F  the  weight, 
NHIST  the  history  number,  and  WC  a  normalization  factor.  I  and  J 12345  are 
indices. 
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J12345  =  1  -  identifies  a  source  particle 

=  2  -  identifies  a  particle  coming  out  of  elastic  scattering 
=  3  -  identifies  a  particle  coming  out  of  inelastic  scattering 
=  4  -  identifies  a  latent  of  the  second  kind  for  FAP  only 
=  5  -  identifies  a  latent  of  the  first  kind  for  FAP  only. 

(In  other  parts  of  the  code  J12345  =  6  identifies  a  transmitted  particle,  J12345  =  7 
an  inelastic  interaction,  and  J 12345  =  8  an  absorption.) 

I  is  irrelevant  (set  to  0)  for  J 12345  =  1,  2,  3  (and  6).  For  FAP  latents  (J  12345  = 
4,  5),  I  is  set  to  IDET,  the  detector  number  for  which  the  latent  applies.  [In  the 
description  of  interactions  (J12345  =  7,  8),  I  is  set  to  IATWT,  a  five-digit  identi¬ 
fier  of  the  element  with  which  the  interaction  occurred.] 

When  the  program  MONTE  calls  the  subroutine  PICK,  it  examines  the  J 12345  ob- 
tained  and  calls  the  proper  routines: 

J 12345  =  1,3  -  both  FAP  and  CARLO 
=  2  -  CARLO  only 
=  4,5  -  FAP  only. 
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9.  PEDIT  -  THE  EDIT  PROGRAM 


The  PEDIT  program  runs  through  the  output  tape,  calculates  average  fluxes  and 
standard  deviations,  and  prints  them  out  in  a  self-explanatory  format.  It  also 
gives  a  tally  of  the  absorptions,  kills,  births,  and  escapes  occurring  in  different 
geometrical  regions. 
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10.  RESPONSE 


The  program  RESPONSE  is  a  separate  program  which  accepts  as  input  the  out¬ 
put  tape  of  an  UNC -SAM-2  problem  and  card  input  to  provide  response  functions 
Rijk(E).  It  computes 

E 

<pik(E)Rij(E)dE 

in  specified  regions  i,  for  all  time  bins  k,  for  all  response  functions  j.  It  also 
computes  associated  statistical  errors.  <pik(E)  is  assumed  to  be  equal  to  the 
average  flux  to  energies  E  in  the  bin. 

An  input  sheet  is  shown  in  Chapter  13. 

.A  first  card  gives: 

NHIST  Number  of  histories  to  be  treated  (must  be  «  the  num¬ 
ber  of  histones  treated  by  the  Monte  Carlo  code). 

The  next  card  gives: 

NRESP  Number  of  response  functions  in  the  problem 
The  description  of  each  response  function  is  as  follows.  A  first  card  gives: 

NEN  Number  of  entries  in  the  table 

NREG  Number  of  regions  in  which  response  is  to  be  calculated 
50H  50  Hollerith  characters  for  identification  purposes. 
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The  next  cards  specify  the  response  function  R(E).  Pairs  of  numbers  E,  E(E) 
ar^*  inputed.  Linear  interpolation  is  assumed  between  entries  in  the  table. 

Finally,  a  list  of  the  (NREG)  re  ;icn  numbers  in  which  the  response  is  to  be  cal 
cuiated  should  be  entered. 


11.  GkSP  -  THE  SECONDARY  PARTICLE  GENERATING  PROGRAM 


This  program  incorporates  simple  interaction  routines  other  than  those  built 
into  the  code.  Given  as  input  points  where  particles  go  into  interaction  written 
on  a  magnetic  tape  during  a  previous  UNC -SAM-2  calculation,  given  a  simplified 
description  of  the  interaction,  and  a  description  of  the  statistical  weights,  GASP 
generates  biased  particles  coming  out  of  collision  and  writes  them  on  another 
magnetic  tape  to  be  used  as  a  source  tape  for  a  future  UNC -SAM-2  calculation. 
The  main  use  of  GASP  is  for  generating  inelastic  and  capture  gamma  rays  due 
to  neution  interaction  in  matter.  It  can  also  be  used  for  generating  the  second 
neutron  from  (n,2n)  reactions,  fission  neutrons,  etc. 

An  input  form  is  shown  in  Chapter  13. 

Input 

A  first  card  specifies: 

NELE  Number  of  elements  present  in  the  configuration,  and  for 
which  GASP  data  are  provided 

IRMAX  Number  of  geometrical  regions  in  the  problem 

IHCUT  Maximum  number  of  histories  to  be  treated 

IGCUT  Maximum  number  of  secondary  particles  to  be  created 

NOUT  Number  of  energy  bins  for  the  edit  of  the  interaction  tape. 


If  NOUT  >  0,  this  is  followed  by  a  specification  of  the  (NOUT  +  1)  energy  bounda- 


ries  of  the  NOUT  energy  bins  for  an  edit  of  the  interaction  tape.  If  NOUT  =  0, 
no  bin  input  is  required  and  the  edit  is  bypassed. 

If  IGCUT  >  0,  this  is  follow  d  by  the  GASP  element  data  for  the  NELE  elements. 

For  each  element,.  this  consists  of: 

1.  Element  Parameter  Card 

IATWT  Integral  part  of  the  atomic  mass  (identification  of  element 
same  as  that  specified  in  DATORG  input). 

la  Number  of  absorption  gamma  output  energies  (0  if 

no  capture  gamma  produced). 

ka  Number  of  corresponding  incident  neutron  energies. 

This  equals  the  number  of  neutron  energy  bins  (0  for 
no  gamma  capture).. 

Lj  Number  of  inelastic  gamma  out  energies. 

Kj  Number  of  corresponding  incident  neutron  energies. 

This  equals  the  number  of  neutron  energy  bins  +  1. 

2.  Capture  Gamma  Energy  Cards  (as  many  as  necessary) 

The  energies  are  given  in  ascending  order  in  eV. 

3.  Capture  Neutron  Energy  Cards  (as  many  as  necessary) 

The  entries  are  in  ascending  order  and  in  eV  and  represent  the  upper 
energy  of  an  energy  bin  (the  low  energy  of  the  lowest  energy  bin  is  not 
entered,  and  is  assumed  to  be  zero). 

4.  Capture  Number  Cards 

Each  entry  represents  the  number  of  gammas  produced  at  each  of  the 
specified  gamma  energies.  These  numbers  are  entered  for  each  neutron 
energy  bin. 

Start  a  new  card  for  each  neutron  bin. 

Repeat  items  2,  3,  and  4  for  inelastic  scattering  input,  if  any. 
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The  element  data  input  is  followed  by  a  description  of  the  region  and  energy  weights 
for  the  secondary  particles,  in  a  format  identical  to  that  of  the  UNO -SAM-2  code. 
Note  that  angular -dependent  weights  are  not  entered. 

If  IGCUT  =  0,  no  element  data,  and  no  weight  specifications  are  required.  The 
generation  of  secondaries  is  bypassed. 

The  following  cards  specify  the  energy  bin  structure  for  the  edit  of  a  source  tape 
(the  source  tape  currently  generated  if  IGCUT  >  0,  or  an  existing  source  tape  of 
IGCUT  =  0). 

A  card  gives: 

NOUT  Number  of  output  bins. 

This  is  followed  by  cards  specifying  the  (NOUT-r-1)  energy  bins  boundaries,  in 
decreasing  order. 

This  is  the  end  of  card  input.  The  code  also  requires  as  input  an  interaction  tape, 
and  a  source  tape  (if  IGCUT^O). 

Output 

The  output  consists  of  the  following. 

1.  If  NOUT  *  0  on  the  first  card,  an  edit  of  the  interaction  tape  is  printed. 
This  lists  the  number  of  interactions  occurring  for  each  geometrical 
region  for  each  element  for  each  energy  bin.  If  NOUT  =  0,  this  edit  is 
bypassed, 

2.  If  IGCUT  *  0,  the  element  data  and  the  weight  specification  are  printed 
back. 
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The  code  then  proceeds  to  generate  secondary  particles  and  writes  them  on  the 
source  tape.  This  is  done  by  reading  a  record  of  a  particle  which  entered  into 
an  interaction.  The  record  consists  of  X,  O,  E,  T,  NHIST,  F,  Wc,  IR,  J 12345, 
IATWT  (where  J12345  is  an  integer  which  specifies  whether  the  interaction  was 
an  inelastic  scattering  or  an  absorption),  and  IATWT  is  a  five-digit  integer  identi¬ 
fier  of  the  element  with  which  the  interaction  occurred.  The  code  proceeds  to 
look  up  the  data  for  element  IATWT  for  the  energies  E y  and  the  number  n(Ey) 
for  each  energy,  of  secondary  particles  to  be  created.  For  each  energy,  it  com¬ 
putes  a  weight 

W  =  WX(IR)  We(E) 
and  computes  the  effective  number 

S(Er)  =n(Ey)£ 

of  secondaries  to  be  picked  from  the  biased  kernel.  The  fractional  part  of  this 
number  is  then  made  either  zero  or  unity  by  generating  a  random  number  £  and 
truncating 

N  -  n  +  £. 

N  directions  are  picked  isotropically,  and  N  source  particles  are  written  on 
the  source  tape,  all  at  the  same  X,  the  same  E,  the  same  T,  the  same  NHIST,  the 
same  F  set  equal  to  the  weight  W,  the  same  Wc,  the  same  IR,  but  with  the  N dif¬ 
ferent  directions  sV. 

When  all  E^’s  have  been  treated  for  the  given  interaction,  the  code  proceeds  to 
the  next,  unless  NHIST  >  IHCUT,  or  the  total  number  of  secondaries  generated 
exceeds  IGCUT,  or  the  last  ("blank”)  interaction  has  been  reached.  At  this  point 
the  code  writes  a  "blank”  source  particle  with  NHIST  =  NHIST  +  1. 
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If  IGCUT  =  0,  all  of  this  calculation  is  bypassed,  but  a  previously  generated  source 
tape  is  assumed  to  be  mounted. 

The  code  then  rewinds  the  source  tape  and  edits  it,  giving  the  number  of  source 
particles  for  each  region  for  each  energy  bin  for  all  source  particles  up  to,  but 
excluding,  the  last  “blank”  one. 

As  can  be  seen,  the  program  GASP  ignores  completely  the  existence  of  possible 
angular-dependent  weights.  If  no  such  biasing  is  used,  GASP  picks  from  the  ex¬ 
act  modified  collision  kernel-  If  angular  importance  is  used  in  the  future  UNC- 
SAM  calculation  for  the  transport  of  the  secondaries,  the  source  particles  are  not 
generated  from  the  completely  modified  kernel.  Hov/ever,  upon  being  picked, 
the  weight  F  =  WxWg  of  the  source  particle  generated  is  divided  by  the  total 
weight  W  =  WxWgWj^  of  the  particle,  giving  an  F  different  from  unity,  and  the 
angular  importance  of  the  particle  picked  will  be  correctly  treated  in  the  trans¬ 
port  problem.  In  other  words,  this  means  that  the  sampling  is  not  as  good  as  it 
could  be,  but  that  the  answers  are  still  correct,  on  the  average. 

The  running  of  the  subsequent  SAM-2  problem  is  straightforward.  The  region 
and  energy  weights  should  be  kept,  if  possible,  the  same  as  those  specified  in 
GASP,  although  the  argument  similar  to  the  one  above  shows  that  there  is  no 
great  harm  in  their  modification. 

The  statistical  aggregate  size  for  the  secondary  problem  should  be  the  same  as, 
or  a  multiple  of,  that  of  the  primary  problem.  This  is  due  to  two  reasons.  If  the 
primary  problem  was  run  using  the  internal  source  generator,  the  energy  of  the 
source  particles  was  picked  in  a  stratified  and  ordered  (decreasing  E)  way. 
Furthermore,  if  the  primary  problem  was  run  with  more  than  one  energy  super¬ 
group,  the  collisions  corresponding  to  a  given  statistical  aggregate  are  not  nec¬ 
essarily  ordered  in  history  number. 


12.  TRANSMIT  -  TREATMENT  OF  TRANSMITTED  PARTICLES 

As  mentioned  earlier  in  the  description  of  CARLO,  one  can  define  a  “transmis¬ 
sion”  region.  Particles  crossing  into  that  region  are  written  on  the  interaction 
tape  with  a  special  index  (J12345=6)  identifying  them. 

The  program  TRANSMIT  is  a  rather  simple  one.  It  runs  through  the  interaction 
tape,  and,  each  time  it  encounters  a  transmitted  particle,  writes  it  out  on  a  source 
tape.  It  also  edits  this  tape  provided  that  output  bins  are  specified.  (See  input 
sheet  in  Chapter  13.) 
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13.  INPUT  FORMS 
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This  Document  Contains 
Missing  Page/s  That  Are 
Unavailable  In  The 
Original  Document 


GiZ,  C'M2j£L- 


HuLmL': 


BEST 

AVAILABLE  COPY 


T 


(CONTINUED) 


GEOM  INPUT  -  (CONTINUED) 

ITYPE  =  1  -  MASQUE  Geometry 

t.  Enter  information  for  I  lines  (1  =  1,  NSRF) 


GEOM  INPUT  -  (CONTINUED) 


This  number  should  be  the  number  of  the  pertinent  region. 


GEOM  INPUT  -  (CONTINUED) 


GEOM  INPUT  -  (CONTINUED) 


Repeat  above  3  cards  (or  each  wedge 


GEOM  INPUT  -  (CONTINUED) 


Repeat  above  3  cards  for  each  non-ordinary  rectangular  parallelepiped 


INPUTD  -  (CONTINUED) 


X 


u 

x 

a 


Si 

u> 


E 

2 


2 

E 


c 

S 

01 

E 


S 

O 


■2 

5 

c 

5 


a 

O 


S«l 371  Ml  3*1401 41 142143  144  1 431441 471 4tl4»|30lftl  1 52  I  S3  I  54 1 53  I  5*1 57 1  3* 1 3* 1*0  I  §11  *2 143 l«4 1*3  I M l«T 


Aiming  Angles  -  enter  NAIML  cards 


RESPONSE  INPUT 


RESPONSE  -  (CONTINUED) 


GASP  INPUT 
General  Specifications 


Enter  Ka  upper  bounds  of  neutron  energy  bins,  in  ascending  order 
(Lower  bound  of  low  energy  bin  is  implied  to  be  0.) 


GASP  -  (CONTINUED) 

Region  Weights  -  enter  NRWL  numbers 


GASP  -  (CONTINUED) 
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